program elgmomc
  use mod_matrix
  use mod_elglike
  use mod_elglike_element
  use mod_elgmc
  use mod_elgdata
  use mod_lutil
  use mod_params
  implicit none

 include 'mpif.h'
 
  character(len=128) :: pfile
  character(len=128) :: datafile, method
  integer :: p1, p2
  integer :: stepsx, stepsy
  type(t_params) :: params
  double precision :: x1, x2, y1, y2


  pfile = "elgmomc.ini"

 
call random_seed()
 
  call MPI_INIT ( global_ierr )
  
  call MPI_COMM_RANK (MPI_COMM_WORLD, global_my_id, global_ierr)
  call MPI_COMM_SIZE (MPI_COMM_WORLD, global_num_procs, global_ierr)
  
  write(*,*) "Hello world! I'm process ", global_my_id, " out of ", global_num_procs, " processes."

  call lutil_get_parameter(21, pfile, "datafile", par_string=datafile)
  call lutil_get_parameter(21, pfile, "lmax", par_int=global_lmax)
  call lutil_get_parameter(21, pfile, "zbins", par_int=zbins)

  call lutil_get_parameter(21, pfile, "method", par_string=method)


  call elgdata_load(global_elgdata, datafile)


  call params_setup(params)

  call lutil_get_parameter(21, pfile, "param_omega_b", par_dp=params%params(1))
  call lutil_get_parameter(21, pfile, "param_omega_nu", par_dp=params%params(2))
  call lutil_get_parameter(21, pfile, "param_omega_w0", par_dp=params%params(3))
  call lutil_get_parameter(21, pfile, "param_omega_wa", par_dp=params%params(4))
  call lutil_get_parameter(21, pfile, "param_omega_A", par_dp=params%params(5))
  

  call lutil_get_parameter(21, pfile, "param_omega_b_steplength", par_dp=params%params_steplength(1))
  call lutil_get_parameter(21, pfile, "param_omega_nu_steplength", par_dp=params%params_steplength(2))
  call lutil_get_parameter(21, pfile, "param_omega_w0_steplength", par_dp=params%params_steplength(3))
  call lutil_get_parameter(21, pfile, "param_omega_wa_steplength", par_dp=params%params_steplength(4))
  call lutil_get_parameter(21, pfile, "param_omega_A_steplength", par_dp=params%params_steplength(5))



  if (global_my_id == global_id_master) then 

!     call lutil_get_parameter(21, pfile, "weightfile", par_string=wfile)

     call lutil_get_parameter(21, pfile, "x1", par_dp=x1)
     call lutil_get_parameter(21, pfile, "x2", par_dp=x2)

     call lutil_get_parameter(21, pfile, "y1", par_dp=y1)
     call lutil_get_parameter(21, pfile, "y2", par_dp=y2)

     call lutil_get_parameter(21, pfile, "stepsx", par_int=stepsx)
     call lutil_get_parameter(21, pfile, "stepsy", par_int=stepsy)

     call lutil_get_parameter(21, pfile, "parameter_1", par_int=p1)
     call lutil_get_parameter(21, pfile, "parameter_2", par_int=p2)



!     call elgmomc_init_master_brute_2D(0.d0, 0.03d0, 0.0d0,0.03d0, global_steps, param_omegaNU, param_wa)
     if (method .eq. "brute_force_2D") then
        call elgmomc_init_master_brute_2D(x1, x2, y1,y2, stepsx, stepsy, p1, p2)
     end if
     if (method .eq. "mcmc") then
        call elgmomc_init_master_mcmc(params)
     end if

  else
     call elgmomc_init_slave()

  endif

  call MPI_FINALIZE ( global_ierr )
  
  stop

end program elgmomc



subroutine elgmomc_init_slave()
  use mod_matrix
  use mod_elglike
  use mod_elglike_element
  implicit none

  include 'mpif.h'
  
  integer :: i = 0
  integer :: ierr, l
  double precision, dimension(10) :: params
  integer, dimension(MPI_STATUS_SIZE*2) :: status

  real*8 :: ll
  real :: likelihood
  type(t_elglike_element) :: elg
  double precision, dimension(2) :: send


  call elglike_element_initialize(elg);



     do while (i == 0)
        call MPI_RECV( params, 10, MPI_DOUBLE_PRECISION, global_id_master, global_tag, MPI_COMM_WORLD, status, global_ierr) 
        l = params(6)
        !write(*,*) "Recieved l " , l
        call elglike_element_setup_matrix(elg, params);
        call elglike_element_calculate_likelihood(elg, params(6), global_elgdata%elements(l));
        likelihood = elg%likelihood
        send(1) = l
        send(2) = likelihood
        call MPI_SEND( send, 2, MPI_DOUBLE_PRECISION, global_id_master, global_tag, MPI_COMM_WORLD,global_ierr) 

    end do


end subroutine elgmomc_init_slave


subroutine elgmomc_init_master_brute_1D()
  use mod_matrix
  use mod_elglike
  use mod_elglike_element
  use mod_params
  implicit none
  type(t_matrix)  :: matrix, invmat
  type(t_elglike) :: elglike
  type(t_params) :: params
  integer :: i,j
  character(len=128) :: filename

  real*8 :: l
  integer :: f

  f = 4
    

  filename = "likelihood.txt"

  params%omegab = 0.2 ! omega m
  params%omeganu = 0.01 ! 0.01 omega_nu
  params%w0 = -1 ! w0
  params%wa = 0.01 ! wa
!  params%amplitude = 1.0 ! amplitude

  call params_setup(params)

  write(*,*) "EH"

  call elglike_initialize(elglike, global_lmax)

  open(unit=f,file=filename,action='write',status='replace')
  do i=1, 20
     write(*,*) "Calculamakting likelihood step ",i
     
     call elglike_calculate_likelihood_master(elglike, params)
     write(f,*) params%params(1), elglike%likelihood
     write(*,*) params%params(1), elglike%likelihood
     params%params(1) = params%params(1) + 0.01
  end do
  close(f)

  write(*,*) "Done!"
  !call MPI_FINALIZE()
  stop
end subroutine elgmomc_init_master_brute_1D

subroutine elgmomc_init_master_brute_2D(X1, X2, Y1, Y2, stepsx, stepsy, p1, p2)
  use mod_matrix
  use mod_elglike
  use mod_elglike_element
  use mod_params
  use mod_likelihood
  implicit none
  integer, intent(in) :: stepsx, stepsy, p1, p2
  double precision, intent(in) :: X1, X2, Y1, Y2
  type(t_matrix)  :: matrix, invmat
  type(t_elglike) :: elglike
  type(t_params) :: params
  integer :: i,j
  character(len=128) :: filename
  double precision :: stepX, stepY
  real*8 :: l
  
  type(t_likelihood) :: likelihood
    

  filename = "likelihood.txt"

  params%omegab = 0.3 ! omega m
  params%omeganu = 0.01 ! 0.01 omega_nu
  params%w0 = -1 ! w0
  params%wa = 0.01 ! wa
!  params%amplitude = 1.0 ! amplitude

  call params_setup(params)

 
  call elglike_initialize(elglike, global_lmax)


  params%params(p1) = x1 ! 0.01 omega_nu
  params%params(p2) = y1 ! w0
  stepX = (x2-x1)/real(stepsx)
  stepY = (y2-y1)/real(stepsy)



  call likelihood_initialize(likelihood, stepsx, stepsy)
  
  do i=1, stepsx

     params%params(p2) = y1 ! w0
     
     do j=1, stepsy
        write(*,*) "Calculamakting likelihood step ",i
        
        call elglike_calculate_likelihood_master(elglike, params)
        !write(f,*) params%params(2), params%params(3), elglike%likelihood
        write(*,*) i,j
        write(*,*) params%params(p1), params%params(p2), elglike%likelihood

        likelihood%val(i,j) = params
        likelihood%val(i,j)%likelihood = elglike%likelihood

        params%params(p2) = params%params(p2) + stepY

     end do
     params%params(p1) = params%params(p1) + stepX
  end do

  call likelihood_save_text(likelihood, filename, p1, p2)


end subroutine elgmomc_init_master_brute_2D


subroutine elgmomc_init_master_mcmc(in_params)
  use mod_matrix
  use mod_elglike
  use mod_elglike_element
  use mod_params
  use mod_elgmc
  implicit none
  type(t_matrix)  :: matrix, invmat
  type(t_elglike) :: elglike
  type(t_params), intent(inout) :: in_params
  type(t_params) :: params
  integer :: i,j
  character(len=128) :: filename

  real*8 :: l

    

  filename = "likelihood.txt"

  

!  call params_setup(params)

  call params_copy(in_params, params)
 
  write(*,*) "Starting mcmc..."
  call elglike_initialize(elglike, global_lmax)

  call elgmc_mcmc(params, 200000, filename) 


  write(*,*) "Done!"
  !call MPI_FINALIZE()
  stop
  
end subroutine elgmomc_init_master_mcmc

  

